#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>

void fft_r2r(double *in,double *out,int size)
{
    fftw_plan_with_nthreads(4);
    fftw_r2r_kind kind = FFTW_REDFT00; // DCT-1
    fftw_plan plan = fftw_plan_r2r_1d(size,in,out,kind,FFTW_ESTIMATE);
    if(plan == NULL) {
        printf("fft:Create fftw_plan error!");
        exit(1);
    }
    fftw_execute(plan);
    fftw_destroy_plan(plan);
    fftw_cleanup_threads();
}

void pfb_ring(double *data,double *subband,double *coeffs,int index,int ntap,int nfft)
{
    double *sum = calloc(sizeof(double),nfft);
    /* index belong to [0,ntap) */
    index %= ntap;
    for(int i = 0;i < nfft;++i){
        for(int j = index;j < ntap;++j){
            sum[i] += data[j * nfft + i] * coeffs[(j-index) * nfft + i];
        }
        for(int j = 0;j < index;++j){
            sum[i] += data[j * nfft + i] * coeffs[(ntap+j-index) * nfft + i];
        }
    }
    for(int i = 0;i < nfft;++i){
        printf("%lf ",sum[i]);
    }
    printf("\n");
    fft_r2r(sum,subband,nfft);
    free(sum);
}

void insert_data(double *in,double *data,int nfft,int ntap,int index)
{
    int iindex = index % ntap;
    for(int i = 0;i < nfft;++i){
        data[iindex * nfft + i] = in[index * nfft + i];
    }
}

int main()
{
    int ntap = 4;
    int nfft = 4;

    double *subband = malloc(sizeof(double) * nfft);

    double *odata = malloc(sizeof(double) * 32);

    double *buffer = calloc(sizeof(double) , 16);

    for(int i = 0;i < 32;++i){
        odata[i] = i+1;
    }

    double coeffs[16] = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};

    insert_data(odata,buffer,nfft,ntap,0);
    insert_data(odata,buffer,nfft,ntap,1);

    // for(int i = 0;i < 16;++i){
    //     printf("%lf ",buffer[i]);
    // }
    printf("\n");
    pfb_ring(odata,subband,coeffs,2,ntap,nfft);

    free(odata);
    free(subband);

    exit(0);
}

